1. Загрузка, очистка по условию

Загрузите датасет very_low_birthweight.RDS (лежит в папке домашнего задания). Это данные о 671 младенце с очень низкой массой тела (<1600 грамм), собранные в Duke University Medical Center доктором Майклом О’Ши c 1981 по 1987 г. Описание переменных см. здесь. Переменными исхода являются колонки ‘dead’, а также время от рождения до смерти или выписки (выводятся из ‘birth’ и ‘exit’. 7 пациентов были выписаны до рождения). Сделайте копию датасета, в которой удалите колонки с количеством пропусков больше 100, а затем удалите все строки с пропусками.

very_low_birthweight <- readRDS("very_low_birthweight.RDS")

full_cols <- very_low_birthweight %>% summarise_all(~ sum(is.na(.))) %>%
                         select (!c(lol, magsulf,meth,toc, pvh, ivh, ipe)) %>%  names () # колонки, где меньше 100 пропусков
very_low_birthweight_с <- very_low_birthweight  %>% select(all_of(full_cols)) %>% 
                                          mutate (id = seq(1,671)) %>%  # введем колонку id для обозначения пациентов до всех удалений
                                                 filter (if_all (where(is.numeric),
                                                                                    ~!is.na(.)))  # удалим строки, где есть NA в численных значениях

# ранее использовала drop.na(), но так учитываюися и категориальные переменные (на 4 кейса меньше остается), а нам все же важнее complete cases по количественным переменным

#summary (very_low_birthweight_с)
#very_low_birthweight_с  %>%  tbl_summary(by = dead) %>% add_p()

very_low_birthweight_с  %>% 
  skimr::skim()
Data summary
Name Piped data
Number of rows 536
Number of columns 20
_______________________
Column type frequency:
factor 4
numeric 16
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
race 3 0.99 FALSE 4 bla: 304, whi: 212, nat: 13, ori: 4
inout 0 1.00 FALSE 2 bor: 449, tra: 87
delivery 1 1.00 FALSE 2 vag: 271, abd: 264
sex 1 1.00 FALSE 2 mal: 268, fem: 267

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
birth 0 1 84.63 1.54 81.51 83.43 84.77 85.82 87.48 ▅▆▇▇▅
exit 0 1 84.76 1.55 81.05 83.56 84.89 85.99 87.72 ▂▆▇▇▅
hospstay 0 1 46.76 63.28 -295.00 21.00 40.00 64.00 797.00 ▁▇▁▁▁
lowph 0 1 7.22 0.13 6.53 7.14 7.22 7.32 7.55 ▁▁▃▇▂
pltct 0 1 204.07 80.70 16.00 146.00 204.00 256.00 571.00 ▃▇▅▁▁
bwt 0 1 1137.02 239.76 400.00 967.50 1160.00 1331.25 1500.00 ▁▃▆▇▇
gest 0 1 29.27 2.23 23.00 28.00 29.00 31.00 38.00 ▂▇▇▁▁
twn 0 1 0.20 0.40 0.00 0.00 0.00 0.00 1.00 ▇▁▁▁▂
apg1 0 1 5.01 2.65 0.00 2.00 5.00 7.00 9.00 ▅▆▅▇▇
vent 0 1 0.54 0.50 0.00 0.00 1.00 1.00 1.00 ▇▁▁▁▇
pneumo 0 1 0.17 0.38 0.00 0.00 0.00 0.00 1.00 ▇▁▁▁▂
pda 0 1 0.20 0.40 0.00 0.00 0.00 0.00 1.00 ▇▁▁▁▂
cld 0 1 0.26 0.44 0.00 0.00 0.00 1.00 1.00 ▇▁▁▁▃
year 0 1 84.63 1.54 81.51 83.44 84.77 85.82 87.48 ▅▆▇▇▅
dead 0 1 0.12 0.33 0.00 0.00 0.00 0.00 1.00 ▇▁▁▁▁
id 0 1 328.72 183.24 2.00 175.50 329.50 479.25 671.00 ▇▇▇▇▆ 
#альтернативный способ расчитать дни в госпитале 
#very_low_birthweight_с %>%  mutate(term = as.numeric(date_decimal(exit)- date_decimal(birth))/86400)

2. Графики плотности распределения, факторы

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

3. Проведите тест на сравнение значений колонки ‘lowph’

# t-test
stat.test <- compare_means(
  lowph ~ inout, data = cleaned_data ,
 method = "wilcox.test"
)
stat.test
## # A tibble: 1 × 8
##   .y.   group1       group2                p     p.adj p.format p.signif method 
##   <chr> <chr>        <chr>             <dbl>     <dbl> <chr>    <chr>    <chr>  
## 1 lowph born at Duke transported 0.000000297 0.0000003 3e-07    ****     Wilcox…
# Create a box plot with rstatix

ggplot(cleaned_data, aes(inout, lowph, color = inout)) + # This is the plot function
  geom_boxplot(width = 0.2, fill = "white", alpha = 0.1) +
  labs(x = "Method", fill = "Method") +
  stat_pvalue_manual(
    stat.test,
    y.position = 7.7, step.increase = 1,
    label = "p.adj"
  )

Н0 - медианы (т.к тест Манн-Уитни) значений lowph для новорожденных, рожденных в Дюке и перевезенных не различаюся. Н0 отвергаем ( p-value очень мало).

Вариант интерпретации - дети транспортированные в госпиталь имеют более низкий показатель lowph (самый низкий pH крови в первые 4 дня жизни) и более низкую выживаемость, чем дети, рожденные в госпитале. Вероятно из других больниц чаще поступают дети в более тяжелом состоянии состоянии, чем рожденные в самом центре. Низкие значения pH крови характерны для тяжелых нарушений гомеостаза, что вероятно требует госпитализации в медицинский центр при университете.

4. Корреляции, графики для корреляций

require(psych)
require (corrplot)

corr_cl_data <- cleaned_data %>% select (where (is.numeric)) # новый датасет в количественными и ранговыми переменными
correlation  <- psych::corr.test(corr_cl_data, method = "spearman", adjust = "fdr") # c поправкой на множ. сравнения

# график матрицы корреляций
corrplot(corr = correlation$r, 
         type="upper", order="hclust", method = "number", col = COL2('PiYG'), 
         p.mat = correlation$p, sig.level = 0.01, insig = "blank")

# network (показывает силу связей и взаимоположение)
library(corrr)
net <- corr_cl_data %>% 
  cor() %>% 
  network_plot()
net

btw и gest значительно скорелированы, hospstay и bwt/gest слабо отрицательно скорелированы.
net plot визуализирует взамосвязи и степень корреляции переменных.

5. Иерархическая кластеризация.

birth_scaled <- scale(corr_cl_data)  # шкалирование 
birth_dist <- dist(birth_scaled)  # вычисление дистанций 

df_dist.hc <- hclust(d = birth_dist,
                     method = "ward.D2")
fviz_dend(df_dist.hc,
          cex = 0.6)
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
##   Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

#library("ape")
#plot(as.phylo(df_dist.hc), type = "fan")

Не очень удачныый способ представления большого количества наблюдений.

6. Одновременный график heatmap и иерархической кластеризации.

library(pheatmap)

hierarchy <- pheatmap(birth_scaled, 
         show_rownames = FALSE, 
         clustering_distance_rows = birth_dist,
         clustering_method = "ward.D2", 
         cutree_rows = 5, #зададим количество кластеров для рядов
         cutree_cols = length(colnames(birth_scaled)),
         angle_col = 45, 
         main = "Dendrograms for clustering rows and columns with heatmap")
hierarchy

На основе дистанций по средним шкалированных значений переменных построили дендрограмму для переменных (сверху). Видно, что дистанция между bwt и gest минимальна (по средним значениям), и эти переменные “похожи” между собой и визуально представляют один кластер, далее скорее всего в один кластер объединены переменные apg1, lowph и pltct. Hostpstay отстоит от остальных переменных максимально далеко по значениям матрицы дистанций. 

По такому же принципу построили дендрограмму для строк, то есть иерархически объединили наблюдения (отдельных новорожденных) по степени близости на матрице дистанций. Значение кол-ва кластеров (отсечки) 5 задали самостоятельно в настройках.

7. PCA анализ

PCA- метод снижения размерности, то есть преобразования данных таким образом, чтобы, снизив количество переменных, сохранить наибольшее число информации об “вариантивности” данных, то есть соблюсти отношения между строками/наблюдениями. Также метод помогает преобразовать данные таким образом, чтобы убрать корреляцию между переменными (в случае “birth” это найденные ранее скоррелированные переменные).

В результате PCA мы получаем “новые” нескоррелрованные переменные - компоненты и их количество равно кол-ву колонок. Компоненты имеют иерархию по признаку того, какой процент дисперсии исходных данных они объясняют (Proportion of variance). Метод позволяет анализировать и строки, и колонки единовременно.

Шкалирование нужно, чтобы корректно соотносить и делать математичекие операции разных переменных друг с другом.

library(factoextra)
library(FactoMineR)
library(ggbiplot) 
## 
## Attaching package: 'ggbiplot'
## The following object is masked from 'package:psych':
## 
##     reflect
birth.pca <- prcomp(corr_cl_data, 
                scale = T) # шкалирование нужно так как датафрейм подаем не шкалированный
summary(birth.pca)
## Importance of components:
##                           PC1    PC2    PC3    PC4     PC5     PC6
## Standard deviation     1.5936 1.0683 0.8825 0.8512 0.73273 0.52813
## Proportion of Variance 0.4233 0.1902 0.1298 0.1208 0.08948 0.04649
## Cumulative Proportion  0.4233 0.6135 0.7433 0.8640 0.95351 1.00000
# birth.pca$rotation

# графическая визуализация доли объясненной изменчивости для каждой из компонент
fviz_eig (birth.pca, addlabels = T, ylim = c (0,50))

В случае датасета birth алгоритм сработал достаточно хорошо, так как 3 компоненты объясняют 74% дисперcии (Cumulative Proportion of variance). 60% для первых 2 компонент означает что в исходных данных много скоррелированных колонок/ либо корреляция сильная.

contr_1 <- fviz_contrib(birth.pca, choice = "var", axes = 1) 
contr_2 <- fviz_contrib(birth.pca, choice = "var", axes = 2) # axes - номер комоненты, из чего состоит каждая компонента?
contr_1

contr_2

pca_graph <- fviz_pca_var(birth.pca, col.var = "contrib")
pca_graph

На графике “Variables- PCA” выше визуализированы 1 и 2 главные компоненты (по осям x и y, % показывает долю объясненной вариации вариативности (дисперсии данных)). Длина стрелки (и цвет стрелки) обозначает сколько корреляции было в переменной и каков вклад переменой в дисперсию, объясненную PC1 и PC2. Видно три группы скоррелированы переменных, причем hospstay и btw/gest отрицательно зависимы. Группа bwt/gest образуют почти прямой угол, что говорит об их вероятной нескоррелированности (по PC2).

На графике “Сontribution” видим из вклада каких переменных состоит компонента. Лучше, когда несколько переменных дают наибольший вклад в объясняемую компонентой дисперсию - в наших данных это наблюдается скорее для PC2 и PC3. Линия на графике показывает значения по оси Y, которое было бы в случае если все переменные вносили одинаковый вклад в компоненту.

8. biplot график для PCA. Цвет по значению колонки ‘dead’.

library(ggbiplot)
biplot <- ggbiplot(birth.pca, 
         scale=0, alpha = 0.3, ellipse = TRUE, groups = cleaned_data$dead) +
         labs(title = "PCA of yield contributing parameters", fill = "dead", color = "dead")
biplot

Высокие показатели для pltct/ apg1/ lowph возможно (направление для поиска) связаны с более благоприятным прогнозом (выживание, dead = 0).

9. 3D vis с помощью plot_ly

    library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
ggplotly(biplot, textposition = "none")
ggplotly(biplot) %>% config(displayModeBar = F)
bipl <- ggbiplot(birth.pca, alpha = 0.2, scale=0, ellipse = TRUE, groups = cleaned_data$dead) +
     geom_point(size = 2,alpha = 0.1, aes(col = cleaned_data$dead, text = paste0("ID: ", cleaned_data$id))) +
     labs(title = "PCA of yield contributing parameters", fill = "dead", color = "dead")
## Warning in geom_point(size = 2, alpha = 0.1, aes(col = cleaned_data$dead, :
## Ignoring unknown aesthetics: text
  labs(title = "PCA of yield contributing parameters", fill = "dead", color = "dead")
## $fill
## [1] "dead"
## 
## $colour
## [1] "dead"
## 
## $title
## [1] "PCA of yield contributing parameters"
## 
## attr(,"class")
## [1] "labels"
ggplotly(bipl, tooltip = "text")

10. Интерпретация PCA анализа

Прописано под графиками в п.7

Почему использовать колонку ‘dead’ для выводов об ассоциации с выживаемостью некорректно?

Так как 2 компоненты , что видно на графике, визуализируют только 2 компоненты, суммарно объясняя 61% дисперсии. Таким образом наблюдаемые тенденции могут служить лишь основанием для выбора дальнейшего направления поиска. Так как анализ эксплораторный, зависимости носят предположительный характер. При высоких значения dim можно бы ожидать что корреляция истинная.

11. UMAP

UMAP позволяет проанализировать близость строк друг другу и плотнее сгруппировать разрозненные наблюдения на основании теории топологии, при снижении размерности. Из-за свойства сохранять локальные расстояния, метод подоходит для анализа отношния наблюдений/строк, но мало подходит для анализа соотношений колонок.

На графике UMAP просележиваются облака точек, однако при стандартных настройках облако исхода 1 (dead) полностью перекрывается с исходом 0, в то время как для PCA перекрытие облаков лишь частичное. PCA дает значительно больше информации, чем UMAP, за счет дополнительной информации о соотношении переменных

C использованием tidymodels:

library(tidymodels)
## ── Attaching packages ────────────────────────────────────── tidymodels 1.2.0 ──
## ✔ broom        1.0.6     ✔ rsample      1.2.1
## ✔ dials        1.3.0     ✔ tune         1.2.1
## ✔ infer        1.0.7     ✔ workflows    1.1.4
## ✔ modeldata    1.4.0     ✔ workflowsets 1.1.0
## ✔ parsnip      1.2.1     ✔ yardstick    1.3.1
## ✔ recipes      1.1.0
## ── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
## ✖ psych::%+%()        masks ggplot2::%+%()
## ✖ scales::alpha()     masks psych::alpha(), ggplot2::alpha()
## ✖ infer::chisq_test() masks rstatix::chisq_test()
## ✖ scales::discard()   masks purrr::discard()
## ✖ plotly::filter()    masks rstatix::filter(), dplyr::filter(), stats::filter()
## ✖ recipes::fixed()    masks stringr::fixed()
## ✖ dials::get_n()      masks rstatix::get_n()
## ✖ dplyr::lag()        masks stats::lag()
## ✖ infer::prop_test()  masks rstatix::prop_test()
## ✖ yardstick::spec()   masks readr::spec()
## ✖ recipes::step()     masks stats::step()
## ✖ infer::t_test()     masks rstatix::t_test()
## • Learn how to get started at https://www.tidymodels.org/start/
library(embed)

umap_p <- recipe (~., data = corr_cl_data) %>% # "тех строка" подаем датасет
    step_normalize(all_predictors()) %>%  # нормализуем
  step_umap(all_predictors()) %>% # проводим umap со станд настройками 15 и 0.01
    prep() %>% # "тех строка" 
    juice() # приводим результаты umap к стандартизированному датасету

umap_p_0 <- umap_p %>% 
    ggplot(aes(UMAP1, UMAP2)) +
    geom_point (aes(color = cleaned_data$dead), 
                            alpha = 0.5, size = 2) +
    labs(color = NULL) +
    theme_minimal()
umap_p_0

Альтернативный способ, пакет (umap)

require (umap)
library (umap)

umap <- umap(corr_cl_data)

df <- data.frame(x = umap$layout[,1],
                 y = umap$layout[,2],
                 dead = cleaned_data$dead)

ggplot(df, aes(x, y, colour = dead)) +
  geom_point()

12. изменение параметров UMAP (n_neighbors и min_dist)

make_umap <- function (corr_cl_data, x,y) { 
  um <- recipe (~., data = corr_cl_data) %>% # "тех строка" подаем датасет
    step_normalize(all_predictors()) %>%  # нормализуем
  step_umap(all_predictors(), neighbors = x, min_dist = y) %>% # проводим umap 
    prep() %>% # "тех строка" 
    juice() # приводим результаты umap к стандартизированному датасету

um_p <- um  %>% 
    ggplot(aes(UMAP1, UMAP2)) +
    geom_point (aes(color = cleaned_data$dead), 
                            alpha = 0.5, size = 2) +
    labs(color = NULL) +
    theme_minimal() +
  geom_text(aes(x = 1, y = 3,
                label = paste0("neigh: ",x,", min_dist: ", y)),
            stat = "unique",
            size = 6, color = "red")
um_p
}
make_umap (corr_cl_data, 3, 0.01)

#make_umap (corr_cl_data, 5, 0.1)
make_umap (corr_cl_data, 50, 0.001)

#make_umap (corr_cl_data, 50, 0.01)

С изменением параметров кол-ва соседей и минмального расстояния изменяется плотность точек, разреженность графика проекции. Самая “плотная” группировка достигается при уменьшении минимальной дистанции и колисечества соседей. Более равномерно распределенная проекуция при увеличении мин. дистанции и кол-ва соседей.

13. пермутация

Сделаем новые колонки на основе данных анализированных ранее.

bwt_permut<- corr_cl_data #clone data for analysis 
bwt_50 <- sample(corr_cl_data$bwt[1:246], 246, replace = FALSE) # сделали вектор на основе первой половину btw без повторов

bwt_permut$bwt_p50 <- c(bwt_50, corr_cl_data$bwt[247:492]) # # пермутировали 50% значений переменной btw без повторов, склеили
bwt_permut$bwt_p100 <- sample(corr_cl_data$bwt, 492, replace = FALSE) # пермутировали 100% значений переменной btw без повторов

PCA для данных с 50% и 100% пермутированными значениями bwt:

# PCA анализ
bwt_50 <- bwt_permut %>% select (!c(bwt, bwt_p100))

bwt_50.pca <- prcomp(bwt_50, 
                scale = T) # шкалирование нужно так как датафрейм подаем не шкалированный
summary(bwt_50.pca)
## Importance of components:
##                          PC1    PC2    PC3    PC4    PC5     PC6
## Standard deviation     1.486 1.0709 0.8876 0.8744 0.8090 0.66171
## Proportion of Variance 0.368 0.1911 0.1313 0.1274 0.1091 0.07298
## Cumulative Proportion  0.368 0.5592 0.6905 0.8179 0.9270 1.00000
biplot_50 <- ggbiplot(bwt_50.pca, 
         scale=0, alpha = 0.3, ellipse = TRUE, groups = cleaned_data$dead) +
         labs(title = "PCA (50% permutated btw)", fill = "dead", color = "dead")

bwt_100 <- bwt_permut %>% select (!c(bwt, bwt_p50))

bwt_100.pca <- prcomp(bwt_100, 
                scale = T) # шкалирование нужно так как датафрейм подаем не шкалированный
summary(bwt_100.pca)
## Importance of components:
##                           PC1    PC2    PC3    PC4    PC5     PC6
## Standard deviation     1.3977 1.0451 1.0064 0.8695 0.8443 0.68746
## Proportion of Variance 0.3256 0.1820 0.1688 0.1260 0.1188 0.07877
## Cumulative Proportion  0.3256 0.5076 0.6764 0.8024 0.9212 1.00000
biplot_100 <- ggbiplot(bwt_100.pca, 
         scale=0, alpha = 0.3, ellipse = TRUE, groups = cleaned_data$dead) +
         labs(title = "PCA (100% permutated btw)", fill = "dead", color = "dead")
biplot_50

biplot_100

ggarrange (biplot, biplot_50 ,biplot_100)

Кулумятивный процент объясненной вариации снижается (на примерно 5 и 10% для 50% и 100% пермутированных переменных соотвественно), что говорит об ухудшении качества модели. Наиболее значительно, ожидаемо, направление стрелки изменилось для самой переменной btw, особенно в случае 100% пермутации. Однако для других переменных, особенно с меньшей вероятностью скореллированных с btw, направление изменилось менее значительно (lowph) или практически не изменилось (hospstay).

Форма облаков и положение точек также претерпели некоторые изменения. По мере увеличения % пермутированнх значений- длина прочих стрелок меняется- они скорее всего теперь берут на себя больший вклад в долю общей обясненной вариации.

UMAP для данных с 50% и 100% пермутированными значениями bwt:

fun_umap <- function (data) {
    umap <- recipe (~., data = data) %>% # "тех строка" подаем датасет
    step_normalize(all_predictors()) %>%  # нормализуем
  step_umap(all_predictors()) %>% # проводим umap со станд настройками
    prep() %>% # "тех строка" 
    juice() # приводим результаты umap к стандартизированному датасету

 umap %>% 
    ggplot(aes(UMAP1, UMAP2)) +
    geom_point (aes(color = cleaned_data$dead), 
                            alpha = 0.5, size = 2) +
    labs(color = NULL) +
    theme_minimal()
}

fun_umap(bwt_50)

fun_umap(bwt_100)

#### 14. Импутация значений и иерархическая кластеризация

# датасет с импутацией значений вместо NA

birth_na_imp<- very_low_birthweight %>% select (hospstay, lowph, pltct, bwt, gest, apg1, dead) %>%
                              mutate_if(is.numeric, ~ (replace(., is.na(.), mean(., na.rm = TRUE)))) %>% 
                              mutate (id = seq(1,671)) %>% relocate (id) %>%
                            select (!c(id, dead)) 
                           


correlation_imp  <- psych::corr.test(birth_na_imp,
                                     method = "spearman", adjust = "fdr") # c поправкой на множ. сравнения

Визуализация взаимосвязей

# график матрицы корреляций
corrplot(corr = correlation_imp$r, 
         type="upper", order="hclust", method = "number", col = COL2('PiYG'), 
         p.mat = correlation_imp$p, sig.level = 0.01, insig = "blank")

# network (показывает силу связей и взаимоположение)
library(corrr)
net_imp <- birth_na_imp %>% 
  cor() %>% 
    network_plot(min_cor = .0)
ggarrange (net, net_imp)

Отмечается то, что пропала большая часть корреляци/взаимосвязи данных. в частности достаточно сильная отрицательная взаимосвязь hospstay (времени в клинике и переменных btw/gest)

Иерархическая кластеризацию на этом датафрейме.

birth_imp_scaled <- scale(birth_na_imp) # шкалирование 
birth_imp_scaled_dist <- dist(birth_imp_scaled)  # вычисление дистанций 

df_dist_imp.hc <- hclust(d = birth_imp_scaled_dist,
                     method = "ward.D2")
fviz_dend(df_dist_imp.hc, horiz = TRUE,
          cex = 0.6)

Одновременный график heatmap и иерархической кластеризации. Интерпретация результата.

library(pheatmap)

hierarchy_imp <- pheatmap(birth_imp_scaled, 
         show_rownames = FALSE, 
         clustering_distance_rows = birth_imp_scaled_dist,
         clustering_method = "ward.D2", 
         cutree_rows = 5, #зададим количество кластеров для рядов
         cutree_cols = length(colnames(birth_imp_scaled)),
         angle_col = 45, 
         main = "Dendrograms for clustering rows and columns with heatmap")

hierarchy_imp

На графике видно как преобразование пропущенных зн методом замены на средние значения очень существенно исказило данные. При том что кластеры по переменным не значимых претерпели изменений, кластеры по рядам изменились существенно. Так же данные потеряли большую долю своей вариативности и усилилась тенденция к среднему, что выражается в большей “одноцветности” хитмепа по сравнению с графиком для датафрейма с удаленными значениями.

Как видно из графиков метод импутации хоть и позволят сохранить нам строки с частично пропущенными значениями, однако существенно искажает взаимосаязи между наблюдениями и может приводить к ложным направлениям поиска. Однако в части случаев, особенно когда данных мало и они невосполнимы, метод замены на среднее дает возможность сохранить информацию в случае неполных кейсов.

15. Импутация значений и PCA/ UMAP

birth.imp.pca <- prcomp(birth_na_imp, 
                scale = T) # шкалирование нужно так как датафрейм подаем не шкалированный
summary(birth.imp.pca)
## Importance of components:
##                           PC1    PC2    PC3    PC4    PC5     PC6
## Standard deviation     1.5368 1.0158 0.9994 0.8468 0.8164 0.47350
## Proportion of Variance 0.3936 0.1720 0.1664 0.1195 0.1111 0.03737
## Cumulative Proportion  0.3936 0.5656 0.7320 0.8516 0.9626 1.00000
# birth.pca$rotation

# графическая визуализация доли объясненной изменчивости для каждой из компонент
fviz_eig (birth.imp.pca, addlabels = T, ylim = c (0,50))

pca_graph_imp <- fviz_pca_var(birth.imp.pca, col.var = "contrib")
pca_graph_imp

biplot_imp <- ggbiplot(birth.imp.pca, 
         scale=0, alpha = 0.3, ellipse = TRUE, groups = very_low_birthweight$dead) +
         labs(title = "PCA (imputed NA)", fill = "dead", color = "dead")
biplot_imp

Процент объясненной дисперсии упал не значительно, для двух первых компонент с 61% до 57%, однако изменился паттерн взаимосвязей, который можно было бы предполагать исходя из графика, полученного в п. 7. Теперь hospstay вносит крайне малый вклад, и нет предполагаемого обратного характера взаимозависимости hospstay и gest/btw.

UMAP

umap_imp <- recipe (~., data = birth_na_imp) %>% # "тех строка" подаем датасет
    step_normalize(all_predictors()) %>%  # нормализуем
  step_umap(all_predictors()) %>% # проводим umap со станд настройками 15 и 0.01
    prep() %>% # "тех строка" 
    juice() # приводим результаты umap к стандартизированному датасету

umap_imp_p<- umap_imp %>% 
    ggplot(aes(UMAP1, UMAP2)) +
    geom_point (aes(color = as.factor(very_low_birthweight$dead)), 
                            alpha = 0.5, size = 2) +
    labs(color = NULL) +
    theme_minimal()
umap_imp_p

Большая часть данных с NA в переменных относилась к умершим детям, поэтому при замене NA на средние точек, соотвествующих погибшим (dead==1) стало значительно больше, и проекция UMAP визуально приобрела большую “полярность” Групп 0 и 1.